% Harmonic balance method and AFT method
 
% W.Q. 2019.08.08

% k1=7;k2=10;c=1;B1=0.3;B2=0.5;gama=0.425;
% Piecewise_Data = [k1;k2;c;B1;B2;gama;counter];
% X = [0;ones(2,1)*1e-5;];
% Relative_Amplitude=HB_Verification(2.8*2*pi,1,Piecewise_Data,X,1)
function   [Relative_Amplitude,Y]=HB_Verification(w,NH,Piecewise_Data,SOLUTIONS_SET,Column_Index)
m=1;                                                                       		% mass parameter
kL=Piecewise_Data(3);                                                                       		% stiffness parameter
cL=Piecewise_Data(4);                                                            % damping parameter
precision=1e-7;                                                            		% Iteration precision

ND=1;
NP=1024;

% Matrices of M, K, C
M=m;
K=kL;
C=cL;                         


% Initialization of X,B,detX
	if((2*NH+1)*ND==size(SOLUTIONS_SET,1))
		X=SOLUTIONS_SET(:,Column_Index);
	else
		X=[SOLUTIONS_SET(:,Column_Index);zeros((2*NH+1)*ND-size(SOLUTIONS_SET,1),1);];    
	end
Y=[X;w];                                                                   		% A hybrid vector is consisted of X and w.

% Initialization of H
   H=ones((2*NH+1)*ND,1)*1e-5;                                             		
   v=1;
   counter=1;
   final_number=50;
%% Correction part start
while (norm(H,2)>precision)
    
	[Hz,~,H,xs,~]=BlackBox2(Piecewise_Data,NH,ND,NP,v,Y,M,C,K);                 % BlackBox2 is called
	detX=-Hz\H;
	X=X+detX;                                                                  	% Update the vector X
	norm(H,2)
	if( counter > final_number)
        break;
    end 
	counter=counter+1;
	Y(1:end-1)=X;                                                             	% Update the vector Y
end
% Correction part  end
    if(norm(H,2)<precision)
        Relative_Amplitude=(max(xs(1:NP))-min(xs(1:NP)))/2;
        Y=[X;w];		
    else
        Relative_Amplitude=0;
    end
end
